#----------------------------------------------------------------------------------------------------#
# Code for Appendixed Analyses
# Engelhardt "Generational Persistence in the Nature of White Racial Attitudes"
#----------------------------------------------------------------------------------------------------#
library(stargazer)
library(lavaan)
library(xtable)
library(plyr)
library(semTools)

# set working directory for folder with data and figures
# setwd("")

# Kam and Burge with Covariates -----------------------------------------------------------
# Loading Data
kb_dat <- read.csv("./data/kam_burge.csv")

# * Traits of Blacks --------------------------------------------------------
RHS <- c("rr_sc*mil + pid7_sc + ideo7_sc + pol_aware + fem + educ_sc + inc_sc")

# Negative
m_negblk <- glm(as.formula(paste0("as.factor(neg_tblk)~", RHS)),
                data = kb_dat, family = binomial(link = "probit"))
summary(m_negblk)

# Positive
m_posblk <- glm(as.formula(paste0("as.factor(pos_tblk)~", RHS)),
                data = kb_dat, family = binomial(link = "probit"))
summary(m_posblk)


# * Individualism -----------------------------------------------------------
RHS <- c("rr_sc*mil + pid7_sc + ideo7_sc + pol_aware + fem + educ_sc + inc_sc")

# Affirmed
m_indvrule <- glm(as.formula(paste0("as.factor(indiv_rule)~", RHS)),
                  data = kb_dat, family = binomial(link = "probit"))
summary(m_indvrule)

# Flouted
m_indvflout <- glm(as.formula(paste0("as.factor(indiv_flouted)~", RHS)),
                   data = kb_dat, family = binomial(link = "probit"))
summary(m_indvflout)

# Broken
m_indvnotengh <- glm(as.formula(paste0("as.factor(indiv_notenough)~", RHS)),
                     data = kb_dat, family = binomial(link = "probit"))
summary(m_indvnotengh)


# * Discrimination --------------------------------------------------------
RHS <- c("rr_sc*mil + pid7_sc + ideo7_sc + pol_aware + fem + educ_sc + inc_sc")

# Exits
m_discexist <- glm(as.formula(paste0("as.factor(discrim_prob_blk)~", RHS)),
                   data = kb_dat, family = binomial(link = "probit"))
summary(m_discexist)

# Denial
m_discnot <- glm(as.formula(paste0("as.factor(discrim_not_prob)~", RHS)),
                 data = kb_dat, family = binomial(link = "probit"))
summary(m_discnot)

# Reverse
m_discrev <- glm(as.formula(paste0("as.factor(discrim_reverse)~", RHS)),
                 data = kb_dat, family = binomial(link = "probit"))
summary(m_discrev)

# * Table -----------------------------------------------------------------
stargazer(m_negblk, m_posblk, m_indvrule, m_indvflout, m_indvnotengh, m_discexist, m_discnot, m_discrev,
          title = "Racial Resentment and Cognitive Responses to Items",
          model.numbers = F, 
          covariate.labels = c("Racial Resentment",
                               "Millennial",
                               "Racial Resentment*Millennial",
                               "Partisanship",
                               "Ideology",
                               "Political Awareness",
                               "Female", "Education", "Income"),
          order = c(1, 2, 9, 3, 4, 5, 6, 7, 8, 10),
          dep.var.caption = "",
          dep.var.labels = c("Negative Traits of Blacks", "Positive Trait of Blacks", 
                             "Individualism Affirmed", "Individualism Flouted", "Individualism Broken",
                             "Discrimination Exists", "Discrimination Denial", "Reverse Discrimination"),
          no.space = T, 
          notes.append = T,
          notes = c("Probit regression results. Standard errors in parentheses."), 
          notes.align = "l",
          out = c("./figures/KB_tab.tex"),
          table.placement = "t",
          digits = 3, omit.stat = c("f", "adj.rsq"), align = T, df = F,
          label = "", header= F
)

# Main Text with Ordered Inputs -------------------------------------------
invar_dat <- subset(anes16, 
                    f2f == 1, 
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")
FIT_invar <- c("chisq.scaled", "df", "cfi.scaled", "srmr", "rmsea.scaled")
VARS <- c("specfavr", "pdisc", "thard", "dless")


GROUP <- "millen"
mod_null <- '
# fixes thresholds across groups
specfavr | c(a1, a1)*t1 + c(a2, a2)*t2 + c(a3, a3)*t3 + c(a4, a4)*t4
pdisc | c(b1, b1)*t1 + c(b2, b2)*t2 + c(b3, b3)*t3 + c(b4, b4)*t4
thard | c(c1, c1)*t1 + c(c2, c2)*t2 + c(c3, c3)*t3 + c(c4, c4)*t4
dless | c(d1, d1)*t1 + c(d2, d2)*t2 + c(d3, d3)*t3 + c(d4, d4)*t4

# fixes variances across groups
specfavr ~ c(t1, t1)*1 
pdisc ~ c(t2, t2)*1 
thard ~ c(t3, t3)*1 
dless ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[FIT_invar]

# Defining invariance model
mod <- '
rr =~ dless + thard + specfavr + pdisc 
# reverse wording
thard ~~ specfavr 
'
fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  ordered = VARS,group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT_invar], 3)
config_fit <- fitMeasures(fit_config)[FIT_invar]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  ordered = VARS, 
                  group = GROUP, 
                  group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT_invar], 3)
metric_fit <- fitMeasures(fit_metric)[FIT_invar]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP, ordered = VARS,
                  group.equal = c("loadings", "thresholds"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT_invar], 3)
scalar_fit <- fitMeasures(fit_scalar)[FIT_invar]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "thresholds", AFIs = fit_measures)
summary(out_scalar)

# * Tables ----------------------------------------------------------------
# Table  
invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          ordered = T,
          title = "Generation Invariance, Face-to-Face Interviews ANES 2016",
          file_name = "./figures/anes_genF2F_models_ordered_tab.tex", 
          varnames = c("Deserve Less", "Try Hard", "Special Favors", "Past Discrimination"))

# Table 
comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit)


permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from face-to-face interviews in the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize_ordered.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/anes_invar_ordered.tex")

# Main Text and Sample Size -----------------------------------------------
invar_dat <- subset(anes16, 
                    f2f == 1, 
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "millen"))

# Defining invariance model
GROUP <- "millen"
FITs <- c("chisq","df", "cfi", "srmr", "rmsea")

replicates <- 1000
config_fit <- array(NA, c(replicates, length(FITs)))
metric_fit <- array(NA, c(replicates, length(FITs)))
scalar_fit <- array(NA, c(replicates, length(FITs)))
colnames(config_fit) <- FITs
colnames(metric_fit) <- FITs
colnames(scalar_fit) <- FITs

mod <- '
rr =~ dless + thard + specfavr + pdisc 
# reverse wording
thard ~~ specfavr 
'

# Bootstrap model
set.seed(1693)
for(i in 1:replicates){
  nmill_dex <- sample(rownames(subset(invar_dat, millen == 0)), 
                      table(invar_dat$millen)[["1"]], replace = F)
  samp_dat <- rbind(subset(invar_dat, millen == 1), invar_dat[nmill_dex,])
  
  fit_config <- cfa(mod, samp_dat, mimic = "Mplus", group = GROUP)
  config_fit[i,] <- fitMeasures(fit_config)[FITs]
  
  fit_metric <- cfa(mod, samp_dat, mimic = "Mplus", group = GROUP, group.equal = c("loadings"))
  metric_fit[i,] <- fitMeasures(fit_metric)[FITs]
  
  fit_scalar <- cfa(mod, samp_dat, mimic = "Mplus", group = GROUP, group.equal = c("loadings", "intercepts"))
  scalar_fit[i,] <- fitMeasures(fit_scalar)[FITs]
  print(paste0((i/replicates)*100,"%"))
}
comb_fit <- rbind(config_fit = apply(config_fit, 2, mean),
                  metric_fit = apply(metric_fit, 2, mean),
                  scalar_fit = apply(scalar_fit, 2, mean))
comb_fit_anes16_f2f <- as.data.frame(comb_fit)
comb_fit_anes16_f2f

pchisq(comb_fit_anes16_f2f$chisq[2]-comb_fit_anes16_f2f$chisq[1],
       comb_fit_anes16_f2f$df[2]-comb_fit_anes16_f2f$df[1], lower.tail = F)
pchisq(comb_fit_anes16_f2f$chisq[3]-comb_fit_anes16_f2f$chisq[2],
       comb_fit_anes16_f2f$df[3]-comb_fit_anes16_f2f$df[2], lower.tail = F)

# Export results for Table
write.csv(comb_fit_anes16_f2f, "./figures/comb_fit_anes16_f2f_sampsize.csv")


# Equivalence in Other Data Collections -----------------------------------

# * 2008 CCAP -------------------------------------------------------------
# Loading Data
ccap08 <- read.csv("./data/ccap08.csv")

# ** Continuous Inputs ----------------------------------------------------
invar_dat <- subset(ccap08, 
                    select = c("M_specfav", "M_deserve", "M_tryhard", "M_pdisc",
                               "millen"))
# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")
VARS <- c("M_specfav", "M_pdisc", "M_tryhard", "M_deserve")

GROUP <- "millen"
mod_null <- '
# fixes intercepts across groups
M_specfav ~~ c(p1, p1)*M_specfav
M_pdisc ~~ c(p2, p2)*M_pdisc
M_tryhard ~~ c(p3, p3)*M_tryhard
M_deserve ~~ c(p4, p4)*M_deserve

# fixes variances across groups
M_specfav ~ c(z1, z1)*1 
M_pdisc ~ c(z2, z2)*1 
M_tryhard ~ c(z3, z3)*1 
M_deserve ~ c(z4, z4)*1 
'

fit_null <- cfa(mod_null, invar_dat, group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
rr =~ M_deserve + M_tryhard + M_specfav + M_pdisc

# reverse wording
M_tryhard ~~ M_specfav
'
fit <- cfa(mod, invar_dat, mimic = "Mplus")
summary(fit, standardized = T)
round(fitMeasures(fit)[FIT], 3)

fit_config <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP, group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric) 

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)


permute_fit <- array(NA, c(3, 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials by Mode",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "ccap08_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{12}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains item intercepts. Data from the 2008 CCAP One residual correlation estimated between \\textit{try hard} and \\textit{special favors}.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/ccap08_invar.tex")

invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          title = "Generation Invariance, 2008 CCAP",
          file_name = "./figures/ccap08_models_tab.tex", 
          varnames = c("Deserve Less", "Try Hard", "Special Favors", "Past Discrimination"))


# ** Ordered Inputs ------------------------------------------------------
invar_dat <- subset(ccap08, 
                    select = c("M_specfav", "M_deserve", "M_tryhard", "M_pdisc",
                               "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")
FIT_invar <- c("chisq.scaled", "df", "cfi.scaled", "srmr", "rmsea.scaled")
VARS <- c("M_specfav", "M_pdisc", "M_tryhard", "M_deserve")

GROUP <- "millen"
mod_null <- '
# fixes thresholds across groups
M_specfav | c(a1, a1)*t1 + c(a2, a2)*t2 + c(a3, a3)*t3 + c(a4, a4)*t4
M_pdisc | c(b1, b1)*t1 + c(b2, b2)*t2 + c(b3, b3)*t3 + c(b4, b4)*t4
M_tryhard | c(c1, c1)*t1 + c(c2, c2)*t2 + c(c3, c3)*t3 + c(c4, c4)*t4
M_deserve | c(d1, d1)*t1 + c(d2, d2)*t2 + c(d3, d3)*t3 + c(d4, d4)*t4

# fixes variances across groups
M_specfav ~ c(z1, z1)*1 
M_pdisc ~ c(z2, z2)*1 
M_tryhard ~ c(z3, z3)*1 
M_deserve ~ c(z4, z4)*1 
'

fit_null <- cfa(mod_null, invar_dat, group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[FIT_invar]

# Defining invariance model 
mod <- '
rr =~ M_deserve + M_tryhard + M_specfav + M_pdisc

# reverse wording
M_tryhard ~~ M_specfav
'
# dless seems best anchor given loading
fit_config <- cfa(mod, invar_dat, mimic = "Mplus", 
                  ordered = VARS, group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT_invar], 3)
config_fit <- fitMeasures(fit_config)[FIT_invar]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  ordered = VARS, group = GROUP, 
                  group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT_invar], 3)
metric_fit <- fitMeasures(fit_metric)[FIT_invar]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  ordered = VARS, group = GROUP, 
                  group.equal = c("loadings", "thresholds"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT_invar], 3)
scalar_fit <- fitMeasures(fit_scalar)[FIT_invar]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "thresholds", AFIs = fit_measures)
summary(out_scalar)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit)

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials by Mode",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_mode_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{12}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains item threshods Data from the 2008 CCAP One residual correlation estimated between \\textit{try hard} and \\textit{special favors}.}}")
print(tab, caption.placement = "top",
      sanitize_ordered.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/ccap08_invar_ordered.tex")

invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          ordered = T,
          title = "Generation Invariance, CCAP 2008",
          file_name = "./figures/ccap08_models_ordered_tab.tex", 
          varnames = c("Deserve Less", "Try Hard", "Special Favors", "Past Discrimination"))






# * 2012 CCAP -------------------------------------------------------------
# Loading Data
ccap12 <- read.csv("./data/ccap12.csv")
# ** Continuous Inputs ----------------------------------------------------
invar_dat <- subset(ccap12, 
                    select = c("pp_deserve", "pp_specfav", "pp_tryhard", "pp_pdisc", "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")


GROUP <- "millen"
mod_null <- '
# fixes intercepts across groups
pp_specfav ~~ c(p1, p1)*pp_specfav
pp_pdisc ~~ c(p2, p2)*pp_pdisc
pp_tryhard ~~ c(p3, p3)*pp_tryhard 
pp_deserve ~~ c(p4, p4)*pp_deserve 

# fixes variances across groups
pp_specfav ~ c(t1, t1)*1 
pp_pdisc ~ c(t2, t2)*1 
pp_tryhard ~ c(t3, t3)*1 
pp_deserve ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
rr =~ pp_deserve + pp_pdisc + pp_specfav + pp_tryhard
#rr ~~ 1*rr
# reverse wording
pp_tryhard ~~ pp_specfav
'
fit <- cfa(mod, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit, standardized = T)
round(fitMeasures(fit)[FIT], 3)


fit_config <- cfa(mod, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
config_fit <- round(fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus", group = GROUP, 
                  group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
metric_fit <- round(fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", group = GROUP, 
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
scalar_fit <- round(fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)
#lavTestScore(fit_metric, epc = TRUE)

## test partial metric equivalence
fit_metric_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("rr =~ pp_specfav"), 
                       group = GROUP, group.equal = c("loadings"))
set.seed(1693) # same permutations
out_metric_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_config, con = fit_metric_part, null = fit_null,
                                 param = "loadings", AFIs = fit_measures)
summary(out_metric_part) 

## test scalar equivalence
fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group.partial = c("rr =~ pp_specfav"),
                  group = GROUP, group.equal = c("loadings", "intercepts"))
scalar_fit <- round(fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric_part, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)
#lavTestScore(fit_scalar, epc = TRUE)

## without freeing special favors
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)
#lavTestScore(fit_scalar, epc = TRUE) 

## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("rr =~ pp_specfav",
                                         "pp_pdisc ~ 1"),
                       group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric_part, con = fit_scalar_part, null = fit_null,
                                 param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part)
#lavTestScore(fit_scalar_part, epc = TRUE) 


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_metric_part)[c("chisq","df", "cfi", "srmr", "rmsea")],
                       fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")])
rownames(comb_fit_part)[4:5] <- c("metric_fit_part",  "scalar_fit_part")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Metric--Partial", "Scalar", "Scalar--Partial")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Metric--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit_part", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric_part@AFI.obs
permute_fit["Metric--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric_part@AFI.pval


permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part@AFI.obs
permute_fit["Scalar--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part@AFI.pval


tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "ccap12_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{13}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from the 2012 CCAP}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/ccap12_invar.tex")


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar,
          metric_partial_1 = fit_metric_part,
          scalar_partial_1 = fit_scalar_part,
          title = "Generation Invariance, 2012 CCAP",
          file_name = "./figures/ccap12_models_tab.tex", 
          varnames = c("Deserve Less", "Past Discrimination", "Try Hard", "Special Favors"))



# ** Ordered Inputs ------------------------------------------------------
invar_dat <- subset(ccap12, 
                    select = c("pp_deserve", "pp_specfav", "pp_tryhard", "pp_pdisc", "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")
FIT_invar <- c("chisq.scaled", "df", "cfi.scaled", "srmr", "rmsea.scaled")
VARS <- c("pp_specfav", "pp_pdisc", "pp_tryhard", "pp_deserve")


GROUP <- "millen"
mod_null <- '
# fixes intercepts across groups
pp_specfav | c(a1, a1)*t1 + c(a2, a2)*t2 + c(a3, a3)*t3 + c(a4, a4)*t4
pp_pdisc | c(b1, b1)*t1 + c(b2, b2)*t2 + c(b3, b3)*t3 + c(b4, b4)*t4
pp_tryhard | c(c1, c1)*t1 + c(c2, c2)*t2 + c(c3, c3)*t3 + c(c4, c4)*t4
pp_deserve | c(d1, d1)*t1 + c(d2, d2)*t2 + c(d3, d3)*t3 + c(d4, d4)*t4

# fixes variances across groups
pp_specfav ~ c(t1, t1)*1 
pp_pdisc ~ c(t2, t2)*1 
pp_tryhard ~ c(t3, t3)*1 
pp_deserve ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[FIT_invar]


# Defining invariance model
mod <- '
rr =~ pp_deserve + pp_pdisc + pp_specfav + pp_tryhard
# reverse wording
pp_tryhard ~~ pp_specfav
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus", 
                  ordered = VARS, group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT_invar], 3)
config_fit <- round(fitMeasures(fit_config)[FIT_invar], 3)

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus", group = GROUP, 
                  ordered = VARS, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT_invar], 3)
metric_fit <- round(fitMeasures(fit_metric)[FIT_invar], 3)

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", group = GROUP, 
                  ordered = VARS, 
                  group.equal = c("loadings", "thresholds"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT_invar], 3)
scalar_fit <- round(fitMeasures(fit_scalar)[FIT_invar], 3)

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "thresholds", AFIs = fit_measures)
summary(out_scalar)

## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("pp_pdisc | t3"),
                       ordered = VARS,
                       group = GROUP, group.equal = c("loadings", "thresholds"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric, con = fit_scalar_part, null = fit_null,
                                 param = "thresholds", AFIs = fit_measures)
summary(out_scalar_part)

## test partial scalar equivalence
fit_scalar_part2 <- cfa(mod, invar_dat, mimic = "Mplus", 
                        group.partial = c("pp_pdisc | t3",
                                          "pp_specfav | t3"),
                        ordered = VARS,
                        group = GROUP, group.equal = c("loadings", "thresholds"))
set.seed(1693) # same permutations
out_scalar_part2 <- permuteMeasEq(nPermute = nperms, 
                                  uncon = fit_metric, con = fit_scalar_part2, null = fit_null,
                                  param = "thresholds", AFIs = fit_measures)
summary(out_scalar_part2)

## test partial scalar equivalence
fit_scalar_part3 <- cfa(mod, invar_dat, mimic = "Mplus", 
                        group.partial = c("pp_pdisc | t3",
                                          "pp_specfav | t3",
                                          "pp_deserve | t3"),
                        ordered = VARS,
                        group = GROUP, group.equal = c("loadings", "thresholds"))
set.seed(1693) # same permutations
out_scalar_part3 <- permuteMeasEq(nPermute = nperms, 
                                  uncon = fit_metric, con = fit_scalar_part3, null = fit_null,
                                  param = "thresholds", AFIs = fit_measures)
summary(out_scalar_part3)

## test partial scalar equivalence
fit_scalar_part4 <- cfa(mod, invar_dat, mimic = "Mplus", 
                        group.partial = c("pp_pdisc | t3",
                                          "pp_specfav | t3",
                                          "pp_deserve | t3",
                                          "pp_tryhard | t4"),
                        ordered = VARS,
                        group = GROUP, group.equal = c("loadings", "thresholds"))
set.seed(1693) # same permutations
out_scalar_part4 <- permuteMeasEq(nPermute = nperms, 
                                  uncon = fit_metric, con = fit_scalar_part4, null = fit_null,
                                  param = "thresholds", AFIs = fit_measures)
summary(out_scalar_part4)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_scalar_part4)[FIT_invar])
rownames(comb_fit_part)[4] <- c("scalar_fit_part4")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar", 
                           "Scalar--Partial4")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial4", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part4", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Scalar--Partial4", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part4@AFI.obs
permute_fit["Scalar--Partial4", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part4@AFI.pval


tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "ccap12_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{13}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from the 2012 CCAP}}")
print(tab, caption.placement = "top",
      sanitize_ordered.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      #comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/ccap12_invar_ordered.tex")

invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          scalar_partial_1 = fit_scalar_part4, 
          ordered = T,
          title = "Generation Invariance, CCAP 2012",
          file_name = "./figures/ccap12_models_ordered_tab.tex", 
          varnames = c("Deserve Less", "Past Discrimination", "Special Favors", "Try Hard"))


# * 2016 Web --------------------------------------------------------------
# ** Continuous Inputs ----------------------------------------------------
invar_dat <- subset(anes16, 
                    f2f == 0,
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "millen"

mod_null <- '
# fixes intercepts across groups
specfavr ~~ c(p1, p1)*specfavr
pdisc ~~ c(p2, p2)*pdisc
thard ~~ c(p3, p3)*thard
dless ~~ c(p4, p4)*dless

# fixes variances across groups
specfavr ~ c(t1, t1)*1 
pdisc ~ c(t2, t2)*1 
thard ~ c(t3, t3)*1 
dless ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
rr =~ pdisc + dless + thard + specfavr 

# reverse wording
thard ~~ specfavr 
'
fit <- cfa(mod, invar_dat, mimic = "Mplus")
summary(fit, standardized = T)
round(fitMeasures(fit)[FIT], 3)

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP,
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)

## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("specfavr ~ 1"),
                       group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric, con = fit_scalar_part, null = fit_null,
                                 param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part)

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")])
rownames(comb_fit_part)[4] <- c("scalar_fit_part")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar", "Scalar--Partial")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part@AFI.obs
permute_fit["Scalar--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part@AFI.pval


tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials vs. Older White Web Respondents",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{past discimrination} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from face-to-face interviews in the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/anes_invar_web.tex")


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar,
          scalar_partial_1 = fit_scalar_part,
          title = "Generation Invariance, Web Interviews ANES 2016",
          file_name = "./figures/anes_web_models_tab.tex", 
          varnames = c("Past Discrimination", "Deserve Less", "Try Hard", "Special Favors"))



# ** Ordered Inputs ------------------------------------------------------
invar_dat <- subset(anes16, 
                    f2f == 0, 
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "millen"))

# Desired Fit Measures
FIT_invar <- c("chisq.scaled", "df", "cfi.scaled", "srmr", "rmsea.scaled")
VARS <- c("specfavr", "pdisc", "thard", "dless")


GROUP <- "millen"
mod_null <- '
# fixes thresholds across groups
specfavr | c(a1, a1)*t1 + c(a2, a2)*t2 + c(a3, a3)*t3 + c(a4, a4)*t4
pdisc | c(b1, b1)*t1 + c(b2, b2)*t2 + c(b3, b3)*t3 + c(b4, b4)*t4
thard | c(c1, c1)*t1 + c(c2, c2)*t2 + c(c3, c3)*t3 + c(c4, c4)*t4
dless | c(d1, d1)*t1 + c(d2, d2)*t2 + c(d3, d3)*t3 + c(d4, d4)*t4

# fixes variances across groups
specfavr ~ c(t1, t1)*1 
pdisc ~ c(t2, t2)*1 
thard ~ c(t3, t3)*1 
dless ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[FIT_invar]


# Defining invariance model
mod <- '
rr =~ dless + thard + specfavr + pdisc 

# reverse wording
thard ~~ specfavr 
'
fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  ordered = VARS, group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT_invar], 3)
config_fit <- fitMeasures(fit_config)[FIT_invar]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  ordered = VARS, group = GROUP, 
                  group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT_invar], 3)
metric_fit <- fitMeasures(fit_metric)[FIT_invar]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  ordered = VARS, group = GROUP,
                  group.equal = c("loadings", "thresholds"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT_invar], 3)
scalar_fit <- fitMeasures(fit_scalar)[FIT_invar]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "thresholds", AFIs = fit_measures)
summary(out_scalar)

## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("specfavr | t3"),
                       ordered = VARS, group = GROUP, 
                       group.equal = c("loadings", "thresholds"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric, con = fit_scalar_part, null = fit_null,
                                 param = "thresholds", AFIs = fit_measures)
summary(out_scalar_part)

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_scalar_part)[FIT_invar])
rownames(comb_fit_part)[4] <- c("scalar_fit_part")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar", "Scalar--Partial")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Scalar--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part@AFI.obs
permute_fit["Scalar--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part@AFI.pval

tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and thresholds to equality. Data from web interviews in the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize_ordered.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/anes_invar_web_ordered.tex")

invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          scalar_partial_1 = fit_scalar_part, 
          ordered = T,
          title = "Generation Invariance, Web Interviews ANES 2016",
          file_name = "./figures/anes_genWeb_models_ordered_tab.tex", 
          varnames = c("Deserve Less", "Try Hard", "Special Favors", "Past Discrimination"))




# * 2018 Pilot ------------------------------------------------------------
# Loading Data
anes18 <- read.csv("./data/anes18.csv")
# ** Continuous Inputs ----------------------------------------------------
invar_dat <- subset(anes18, 
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "millen"
mod_null <- '
# fixes intercepts across groups
specfavr ~~ c(p1, p1)*specfavr
pdisc ~~ c(p2, p2)*pdisc
thard ~~ c(p3, p3)*thard
dless ~~ c(p4, p4)*dless

# fixes variances across groups
specfavr ~ c(t1, t1)*1 
pdisc ~ c(t2, t2)*1 
thard ~ c(t3, t3)*1 
dless ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]

# Defining invariance model
mod <- '
rr =~ pdisc + dless + thard + specfavr

# reverse wording
thard ~~ specfavr 
'
fit <- cfa(mod, invar_dat, mimic = "Mplus")
summary(fit, standardized = T)
round(fitMeasures(fit)[FIT], 3)

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP,
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)

permute_fit <- array(NA, c(nrow(comb_fit), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{13}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{past discrimination} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from the 2018 ANES pilot.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      #comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/anes18_invar.tex")


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar,
          title = "Generation Invariance, ANES 2018 Pilot",
          file_name = "./figures/anes_18_models_tab.tex", 
          varnames = c("Past Discrimination", "Deserve Less", "Try Hard", "Special Favors"))


# ** Ordered Inputs ------------------------------------------------------
invar_dat <- subset(anes18, 
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "millen"))

# Desired Fit Measures
FIT_invar <- c("chisq.scaled", "df", "cfi.scaled", "srmr", "rmsea.scaled")
VARS <- c("specfavr", "pdisc", "thard", "dless")

GROUP <- "millen"
mod_null <- '
# fixes thresholds across groups
specfavr | c(a1, a1)*t1 + c(a2, a2)*t2 + c(a3, a3)*t3 + c(a4, a4)*t4
pdisc | c(b1, b1)*t1 + c(b2, b2)*t2 + c(b3, b3)*t3 + c(b4, b4)*t4
thard | c(c1, c1)*t1 + c(c2, c2)*t2 + c(c3, c3)*t3 + c(c4, c4)*t4
dless | c(d1, d1)*t1 + c(d2, d2)*t2 + c(d3, d3)*t3 + c(d4, d4)*t4

# fixes variances across groups
specfavr ~ c(t1, t1)*1 
pdisc ~ c(t2, t2)*1 
thard ~ c(t3, t3)*1 
dless ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[FIT_invar]

# Defining invariance model
mod <- '
rr =~ dless + thard + specfavr + pdisc 

# reverse wording
thard ~~ specfavr 
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, ordered = VARS)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT_invar], 3)
config_fit <- fitMeasures(fit_config)[FIT_invar]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, ordered = VARS, 
                  group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT_invar], 3)
metric_fit <- fitMeasures(fit_metric)[FIT_invar]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP, ordered = VARS,
                  group.equal = c("loadings", "thresholds"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT_invar], 3)
scalar_fit <- fitMeasures(fit_scalar)[FIT_invar]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "thresholds", AFIs = fit_measures)
summary(out_scalar)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- comb_fit

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{13}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and thresholds to equality. Data from the 2018 ANES pilot.}}")
print(tab, caption.placement = "top",
      sanitize_ordered.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      #comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/anes18_invar_ordered.tex")

invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar,
          ordered = T,
          title = "Generation Invariance, 2018 ANES Pilot",
          file_name = "./figures/anes18_models_ordered_tab.tex", 
          varnames = c("Deserve Less", "Try Hard", "Special Favors","Past Discrimination"))



# Descriptive Differences in Kam and Burge --------------------------------
# Loading Data
kb_dat <- read.csv("./data/kam_burge.csv")

# * Traits of Blacks --------------------------------
# Percentages 
prop.table(table(kb_dat$mil, kb_dat$neg_tblk), 1)
# hypothesis test
chisq.test(table(kb_dat$mil, kb_dat$neg_tblk))

# Percentages
prop.table(table(kb_dat$mil, kb_dat$pos_tblk), 1)
# hypothesis test
chisq.test(table(kb_dat$mil, kb_dat$pos_tblk))

# * Individualism --------------------------------
# Percentages
prop.table(table(kb_dat$mil, kb_dat$indiv_rule), 1)
# hypothesis test
chisq.test(table(kb_dat$mil, kb_dat$indiv_rule))
# Percentages
prop.table(table(kb_dat$mil, kb_dat$indiv_notenough), 1)
# hypothesis test
chisq.test(table(kb_dat$mil, kb_dat$indiv_notenough))
# Percentages
prop.table(table(kb_dat$mil, kb_dat$indiv_flouted), 1)
# hypothesis test
chisq.test(table(kb_dat$mil, kb_dat$indiv_flouted))

# * Discrimination --------------------------------
# Percentages
prop.table(table(kb_dat$mil, kb_dat$discrim_prob_blk), 1)
# hypothesis test
chisq.test(table(kb_dat$mil, kb_dat$discrim_prob_blk))
# Percentages
prop.table(table(kb_dat$mil, kb_dat$discrim_not_prob), 1)
# hypothesis test
chisq.test(table(kb_dat$mil, kb_dat$discrim_not_prob))
# Percentages
prop.table(table(kb_dat$mil, kb_dat$discrim_reverse), 1)
# hypothesis test
chisq.test(table(kb_dat$mil, kb_dat$discrim_reverse))

# Mode Comparison for Millennials -----------------------------------------
# * Continuous Inputs -----------------------------------------------------
invar_dat <- subset(anes16, 
                    millen == 1, 
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "f2f"))
# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")
VARS <- c("specfavr", "pdisc", "thard", "dless")

GROUP <- "f2f"
mod_null <- '
# fixes intercepts across groups
specfavr ~~ c(p1, p1)*specfavr
pdisc ~~ c(p2, p2)*pdisc
thard ~~ c(p3, p3)*thard
dless ~~ c(p4, p4)*dless

# fixes variances across groups
specfavr ~ c(z1, z1)*1 
pdisc ~ c(z2, z2)*1 
thard ~ c(z3, z3)*1 
dless ~ c(z4, z4)*1 
'

fit_null <- cfa(mod_null, invar_dat, group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
rr =~ dless + thard + specfavr + pdisc
# reverse wording
thard ~~ specfavr 
'
fit <- cfa(mod, invar_dat, mimic = "Mplus")
summary(fit, standardized = T)
round(fitMeasures(fit)[FIT], 3)

fit_config <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP, group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar,
          title = "Mode Invariance, Millennial Whites ANES 2016",
          file_name = "./figures/anes_millen_models_tab.tex", 
          varnames = c("Deserve Less", "Try Hard", "Special Favors", "Past Discrimination"))


permute_fit <- array(NA, c(3, 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials by Mode",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_mode_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{12}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains item intercepts. Data from the 2016 ANES. One residual correlation estimated between \\textit{try hard} and \\textit{special favors}.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/anes_mode_invar.tex")


# * Ordered Inputs --------------------------------------------------------
invar_dat <- subset(anes16, 
                    millen == 1,
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "f2f"))
# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")
FIT_invar <- c("chisq.scaled", "df", "cfi.scaled", "srmr", "rmsea.scaled")
VARS <- c("specfavr", "pdisc", "thard", "dless")

GROUP <- "f2f"
mod_null <- '
# fixes thresholds across groups
specfavr | c(a1, a1)*t1 + c(a2, a2)*t2 + c(a3, a3)*t3 + c(a4, a4)*t4
pdisc | c(b1, b1)*t1 + c(b2, b2)*t2 + c(b3, b3)*t3 + c(b4, b4)*t4
thard | c(c1, c1)*t1 + c(c2, c2)*t2 + c(c3, c3)*t3 + c(c4, c4)*t4
dless | c(d1, d1)*t1 + c(d2, d2)*t2 + c(d3, d3)*t3 + c(d4, d4)*t4

# fixes variances across groups
specfavr ~ c(z1, z1)*1 
pdisc ~ c(z2, z2)*1 
thard ~ c(z3, z3)*1 
dless ~ c(z4, z4)*1 
'

fit_null <- cfa(mod_null, invar_dat, group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[FIT_invar]


# Defining invariance model
mod <- '
rr =~ dless + thard + specfavr + pdisc

# rr ~~ 1*rr

# reverse wording
thard ~~ specfavr # needed?
'
fit_config <- cfa(mod, invar_dat, mimic = "Mplus", 
                  ordered = VARS, group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT_invar], 3)
config_fit <- fitMeasures(fit_config)[FIT_invar]
# modindices(fit_config)

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  ordered = VARS, group = GROUP, 
                  group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT_invar], 3)
metric_fit <- fitMeasures(fit_metric)[FIT_invar]
# modindices(fit_metric)

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  ordered = VARS, group = GROUP, 
                  group.equal = c("loadings", "thresholds"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT_invar], 3)
scalar_fit <- fitMeasures(fit_scalar)[FIT_invar]
# modindices(fit_scalar)

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "thresholds", AFIs = fit_measures)
summary(out_scalar)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit)

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq.scaled", "cfi.scaled", "srmr", "rmsea.scaled")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval


tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials by Mode",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_mode_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{12}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains item intercepts. Data from the 2016 ANES. One residual correlation estimated between \\textit{try hard} and \\textit{special favors}.}}")
print(tab, caption.placement = "top",
      sanitize_ordered.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/anes_mode_invar_ordered.tex")



invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          ordered = T,
          title = "Mode Invariance, Millennials ANES 2016",
          file_name = "./figures/anes_mode_models_ordered_tab.tex", 
          varnames = c("Deserve Less", "Try Hard", "Special Favors", "Past Discrimination"))

# Additional Construct Equivalence ----------------------------------------

# * Moral Traditionalism --------------------------------------------------
# ** F2F -------------------------------------------------------------------
invar_dat <- subset(anes16, 
                    f2f == 1,
                    select = c("world_change", "lifestyles", "tol_other_morals", "emph_fam", 
                               "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "millen"
mod_null <- '
# fixes intercepts across groups
world_change ~~ c(p1, p1)*world_change
lifestyles ~~ c(p2, p2)*lifestyles
tol_other_morals ~~ c(p3, p3)*tol_other_morals
emph_fam ~~ c(p4, p4)*emph_fam

# fixes variances across groups
world_change ~ c(t1, t1)*1 
lifestyles ~ c(t2, t2)*1 
tol_other_morals ~ c(t3, t3)*1 
emph_fam ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
mtrad =~ tol_other_morals + lifestyles + world_change + emph_fam 

# reverse wording
lifestyles ~~ emph_fam 
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP,
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            #null = fit_null, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)


## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)

## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("emph_fam ~ 1"),
                       group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric, con = fit_scalar_part, null = fit_null,
                                 param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")])
rownames(comb_fit_part)[4] <- c("scalar_fit_part")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar", "Scalar--Partial")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part@AFI.obs
permute_fit["Scalar--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part@AFI.pval


tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Moral Traditionalism, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{tolerate other morals} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from face-to-face interviews in the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      #comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/mtrad_anes_invar_f2f.tex")


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          scalar_partial_1 = fit_scalar_part, 
          title = "Generation Invariance, Face-to-Face ANES 2016",
          file_name = "./figures/mtrad_anes_f2f_models_tab.tex", 
          varnames = c("Tolerate Others", "Lifestyles", "World changing", "Family Ties"))

# *** Substantive Differences ----------------------------------------------
MODEL <- fit_config
pars <- parameterEstimates(MODEL)
SDI2.UDI2(4, 
          pars$est[which(pars$op == "=~" & pars$group == 1)],
          pars$est[which(pars$op == "=~" & pars$group == 2)],
          pars$est[which(pars$op == "~1" & pars$lhs != "mtrad" & pars$group == 1)],
          pars$est[which(pars$op == "~1" & pars$lhs != "mtrad" & pars$group == 2)],
          sqrt(MODEL@Data@ov$var),
          pars$est[which(pars$op == "~1" & pars$lhs == "mtrad" & pars$group == 2)],
          sqrt(pars$est[which(pars$op == "~~" & pars$lhs == "mtrad" & pars$group == 2)]))

# ** Web -------------------------------------------------------------------
invar_dat <- subset(anes16, 
                    f2f == 0,
                    select = c("world_change", "lifestyles", "tol_other_morals", "emph_fam", 
                               "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "millen"
mod_null <- '
# fixes intercepts across groups
world_change ~~ c(p1, p1)*world_change
lifestyles ~~ c(p2, p2)*lifestyles
tol_other_morals ~~ c(p3, p3)*tol_other_morals
emph_fam ~~ c(p4, p4)*emph_fam

# fixes variances across groups
world_change ~ c(t1, t1)*1 
lifestyles ~ c(t2, t2)*1 
tol_other_morals ~ c(t3, t3)*1 
emph_fam ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
mtrad =~ lifestyles + tol_other_morals + world_change + emph_fam 

# reverse wording
lifestyles ~~ emph_fam 
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP,
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)


## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)

## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("emph_fam ~ 1"),
                       group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric, con = fit_scalar_part, null = fit_null,
                                 param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part)

## test partial scalar equivalence
fit_scalar_part2 <- cfa(mod, invar_dat, mimic = "Mplus", 
                        group.partial = c("emph_fam ~ 1",
                                          "tol_other_morals ~ 1"),
                        group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part2 <- permuteMeasEq(nPermute = nperms, 
                                  uncon = fit_metric, con = fit_scalar_part2, null = fit_null,
                                  param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part2)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")],
                       fitMeasures(fit_scalar_part2)[c("chisq","df", "cfi", "srmr", "rmsea")])
rownames(comb_fit_part)[4:5] <- c("scalar_fit_part", "scalar_fit_part2")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar", "Scalar--Partial", "Scalar--Partial2")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part@AFI.obs
permute_fit["Scalar--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part@AFI.pval

permute_fit["Scalar--Partial2", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part2", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial2", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part2@AFI.obs
permute_fit["Scalar--Partial2", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part2@AFI.pval


tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Moral Traditionalism, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{lifestyles} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from web interviews in the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      #comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/mtrad_anes_invar_web.tex")


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          scalar_partial_1 = fit_scalar_part, 
          scalar_partial_2 = fit_scalar_part2, 
          title = "Generation Invariance, Web ANES 2016",
          file_name = "./figures/mtrad_anes_web_models_tab.tex", 
          varnames = c("Lifestyles", "Tolerate Others", "World changing", "Family Ties"))

# *** Substantive Differences ----------------------------------------------
MODEL <- fit_config
pars <- parameterEstimates(MODEL)
SDI2.UDI2(4, 
          pars$est[which(pars$op == "=~" & pars$group == 1)],
          pars$est[which(pars$op == "=~" & pars$group == 2)],
          pars$est[which(pars$op == "~1" & pars$lhs != "mtrad" & pars$group == 1)],
          pars$est[which(pars$op == "~1" & pars$lhs != "mtrad" & pars$group == 2)],
          sqrt(MODEL@Data@ov$var),
          pars$est[which(pars$op == "~1" & pars$lhs == "mtrad" & pars$group == 2)],
          sqrt(pars$est[which(pars$op == "~~" & pars$lhs == "mtrad" & pars$group == 2)]))

# ** Mode --------------------------------------------------------------------
invar_dat <- subset(anes16, 
                    millen == 1,
                    select = c("world_change", "lifestyles", "tol_other_morals", "emph_fam", 
                               "f2f"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "f2f"
mod_null <- '
# fixes intercepts across groups
world_change ~~ c(p1, p1)*world_change
lifestyles ~~ c(p2, p2)*lifestyles
tol_other_morals ~~ c(p3, p3)*tol_other_morals
emph_fam ~~ c(p4, p4)*emph_fam

# fixes variances across groups
world_change ~ c(t1, t1)*1 
lifestyles ~ c(t2, t2)*1 
tol_other_morals ~ c(t3, t3)*1 
emph_fam ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
mtrad =~ emph_fam + world_change + tol_other_morals + lifestyles 

# reverse wording
lifestyles ~~ emph_fam 
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP,
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)


## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)


## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("tol_other_morals ~ 1"),
                       group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric, con = fit_scalar_part, null = fit_null,
                                 param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")])
rownames(comb_fit_part)[4] <- c("scalar_fit_part")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar", "Scalar--Partial")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part@AFI.obs
permute_fit["Scalar--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part@AFI.pval

tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Equivalence of Moral Traditionalism, Mode",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{emphasize family} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      #comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/mtrad_anes_invar_mode.tex")


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          scalar_partial_1 = fit_scalar_part, 
          title = "Mode Equivalence,ANES 2016",
          file_name = "./figures/mtrad_anes_mode_models_tab.tex", 
          varnames = c("Family Ties", "World changing", "Tolerate Others", "Lifestyles"))

# *** Substantive Differences ----------------------------------------------
MODEL <- fit_config
pars <- parameterEstimates(MODEL)
SDI2.UDI2(4, 
          pars$est[which(pars$op == "=~" & pars$group == 1)],
          pars$est[which(pars$op == "=~" & pars$group == 2)],
          pars$est[which(pars$op == "~1" & pars$lhs != "mtrad" & pars$group == 1)],
          pars$est[which(pars$op == "~1" & pars$lhs != "mtrad" & pars$group == 2)],
          sqrt(MODEL@Data@ov$var),
          pars$est[which(pars$op == "~1" & pars$lhs == "mtrad" & pars$group == 2)],
          sqrt(pars$est[which(pars$op == "~~" & pars$lhs == "mtrad" & pars$group == 2)]))


# * Egalitarianism --------------------------------------------------------
# ** F2F -------------------------------------------------------------------
invar_dat <- subset(anes16, 
                    f2f == 1,
                    select = c("ign_eql_opp", "eql_opp", "more_chance", "treat_fair", 
                               "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "millen"
mod_null <- '
# fixes intercepts across groups
ign_eql_opp ~~ c(p1, p1)*ign_eql_opp
eql_opp ~~ c(p2, p2)*eql_opp
more_chance ~~ c(p3, p3)*more_chance
treat_fair ~~ c(p4, p4)*treat_fair

# fixes variances across groups
ign_eql_opp ~ c(t1, t1)*1 
eql_opp ~ c(t2, t2)*1 
more_chance ~ c(t3, t3)*1 
treat_fair ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
egal =~ ign_eql_opp + more_chance + eql_opp + treat_fair 

# reverse wording
eql_opp ~~ treat_fair
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP,
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)

## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("eql_opp ~ 1"),
                       group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric, con = fit_scalar_part, null = fit_null,
                                 param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")])
rownames(comb_fit_part)[4] <- c("scalar_fit_part")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar", "Scalar--Partial")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part@AFI.obs
permute_fit["Scalar--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part@AFI.pval


tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Egalitarianism, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{ignore equal opportunity} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from face-to-face interviews in the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      #comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/egal_anes_invar_f2f.tex")


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          scalar_partial_1 = fit_scalar_part, 
          title = "Generation Invariance, Face-to-Face ANES 2016",
          file_name = "./figures/egal_anes_f2f_models_tab.tex", 
          varnames = c("ign_eql_opp", "more_chance", "eql_opp", "treat_fair"))

# *** Substantive Differences ----------------------------------------------
MODEL <- fit_config
pars <- parameterEstimates(MODEL)
SDI2.UDI2(4, 
          pars$est[which(pars$op == "=~" & pars$group == 1)],
          pars$est[which(pars$op == "=~" & pars$group == 2)],
          pars$est[which(pars$op == "~1" & pars$lhs != "egal" & pars$group == 1)],
          pars$est[which(pars$op == "~1" & pars$lhs != "egal" & pars$group == 2)],
          sqrt(MODEL@Data@ov$var),
          pars$est[which(pars$op == "~1" & pars$lhs == "egal" & pars$group == 2)],
          sqrt(pars$est[which(pars$op == "~~" & pars$lhs == "egal" & pars$group == 2)]))

# ** Web -------------------------------------------------------------------
invar_dat <- subset(anes16, 
                    f2f == 0,
                    select = c("ign_eql_opp", "eql_opp", "more_chance", "treat_fair", 
                               "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "millen"
mod_null <- '
# fixes intercepts across groups
ign_eql_opp ~~ c(p1, p1)*ign_eql_opp
eql_opp ~~ c(p2, p2)*eql_opp
more_chance ~~ c(p3, p3)*more_chance
treat_fair ~~ c(p4, p4)*treat_fair

# fixes variances across groups
ign_eql_opp ~ c(t1, t1)*1 
eql_opp ~ c(t2, t2)*1 
more_chance ~ c(t3, t3)*1 
treat_fair ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
egal =~ ign_eql_opp + more_chance + eql_opp + treat_fair 

# reverse wording
eql_opp ~~ treat_fair 
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP,
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            #null = fit_null, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)


## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)

## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("eql_opp ~ 1"),
                       group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric, con = fit_scalar_part, null = fit_null,
                                 param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part)

## test partial scalar equivalence
fit_scalar_part2 <- cfa(mod, invar_dat, mimic = "Mplus", 
                        group.partial = c("eql_opp ~ 1",
                                          "treat_fair ~ 1"),
                        group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part2 <- permuteMeasEq(nPermute = nperms, 
                                  uncon = fit_metric, con = fit_scalar_part2, null = fit_null,
                                  param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part2)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")],
                       fitMeasures(fit_scalar_part2)[c("chisq","df", "cfi", "srmr", "rmsea")])
rownames(comb_fit_part)[4:5] <- c("scalar_fit_part", "scalar_fit_part2")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar", "Scalar--Partial", "Scalar--Partial2")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part@AFI.obs
permute_fit["Scalar--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part@AFI.pval

permute_fit["Scalar--Partial2", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part2", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial2", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part2@AFI.obs
permute_fit["Scalar--Partial2", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part2@AFI.pval


tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Egalitarianism, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{ignore qual opportunity} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from web interviews in the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      #comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/egal_anes_invar_web.tex")


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          scalar_partial_1 = fit_scalar_part, 
          scalar_partial_2 = fit_scalar_part2, 
          title = "Generation Invariance, Web ANES 2016",
          file_name = "./figures/egal_anes_web_models_tab.tex", 
          varnames = c("ign_eql_opp","more_chance", "eql_opp", "treat_fair"))

# *** Substantive Differences ----------------------------------------------
MODEL <- fit_config
pars <- parameterEstimates(MODEL)
SDI2.UDI2(4, 
          pars$est[which(pars$op == "=~" & pars$group == 1)],
          pars$est[which(pars$op == "=~" & pars$group == 2)],
          pars$est[which(pars$op == "~1" & pars$lhs != "egal" & pars$group == 1)],
          pars$est[which(pars$op == "~1" & pars$lhs != "egal" & pars$group == 2)],
          sqrt(MODEL@Data@ov$var),
          pars$est[which(pars$op == "~1" & pars$lhs == "egal" & pars$group == 2)],
          sqrt(pars$est[which(pars$op == "~~" & pars$lhs == "egal" & pars$group == 2)]))


# * Mode --------------------------------------------------------------------
invar_dat <- subset(anes16, 
                    millen == 1,
                    select = c("ign_eql_opp", "eql_opp", "more_chance", "treat_fair", 
                               "f2f"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "f2f"
mod_null <- '
# fixes intercepts across groups
ign_eql_opp ~~ c(p1, p1)*ign_eql_opp
eql_opp ~~ c(p2, p2)*eql_opp
more_chance ~~ c(p3, p3)*more_chance
treat_fair ~~ c(p4, p4)*treat_fair

# fixes variances across groups
ign_eql_opp ~ c(t1, t1)*1 
eql_opp ~ c(t2, t2)*1 
more_chance ~ c(t3, t3)*1 
treat_fair ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
egal =~ ign_eql_opp + more_chance + eql_opp + treat_fair 

# reverse wording
eql_opp ~~ treat_fair 
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP,
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)


## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)


## test partial scalar equivalence
fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("treat_fair ~ 1"),
                       group = GROUP, group.equal = c("loadings", "intercepts"))
set.seed(1693) # same permutations
out_scalar_part <- permuteMeasEq(nPermute = nperms, 
                                 uncon = fit_metric, con = fit_scalar_part, null = fit_null,
                                 param = "intercepts", AFIs = fit_measures)
summary(out_scalar_part)


comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_part <- rbind(comb_fit,
                       fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")])
rownames(comb_fit_part)[4] <- c("scalar_fit_part")

permute_fit <- array(NA, c(nrow(comb_fit_part), 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar", "Scalar--Partial")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["config_fit", c("chisq", "cfi", "srmr", "rmsea")]

permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

permute_fit["Scalar--Partial", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit_part["scalar_fit_part", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar--Partial", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar_part@AFI.obs
permute_fit["Scalar--Partial", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar_part@AFI.pval


tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Egalitarianism, Mode",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{ignore equal opportunity} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      #comment = F,
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/egal_anes_invar_mode.tex")


invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          scalar_partial_1 = fit_scalar_part, 
          title = "Mode Invariance,ANES 2016",
          file_name = "./figures/egal_anes_mode_models_tab.tex", 
          varnames = c("ign_eql_opp", "more_chance", "eql_opp", "treat_fair"))


